#include "feMeshFunction.h"
#include "feMesh.h"
#include <algorithm>
#include <math.h>
#include "feInterpolation.h"
#include "feObserver.h"
#include <fstream> 



feSharedPtr<feMeshFunction> feSimple1dMeshFunction::DiffDtOperation::operator() (feMeshFunction *iFunction)
{
	//cast
	feSimple1dMeshFunction *f = dynamic_cast<feSimple1dMeshFunction*>(iFunction);
	if(!f) return 0;


	//allocate new function
	feSharedPtr<feMeshFunction> df(new feSimple1dMeshFunction(f->Mesh));
	df->Clear();

	
	//fill new function
	feVector diflayer(f->Mesh->NodesCount());
	for(int i=0; i<f->LayersCount(); i++)
	{
		//calc du/dt
		for(int j=0; j<f->Mesh->NodesCount(); j++)		
			diflayer[j] = f->GetDuDt((int)j, (int)i);
		
		//add layer
		df->AddLayer(f->GetTime(i), &diflayer);
	}

	return df;
}


//---------------------------------------------------------------------------------------------------------------------------------
feSharedPtr<feMeshFunction> feSimple1dMeshFunction::SubOperation::operator() (feMeshFunction *iFunction1, feMeshFunction *iFunction2)
{
	//cast
	feSimple1dMeshFunction *f1 = dynamic_cast<feSimple1dMeshFunction*>(iFunction1);
	feSimple1dMeshFunction *f2 = dynamic_cast<feSimple1dMeshFunction*>(iFunction2);
	
	if(!f1 || !f2) return 0;
	if(f1->GetMesh()->NodesCount() != f2->GetMesh()->NodesCount()) return 0;
	if(f1->LayersCount() != f2->LayersCount()) return 0;


	//allocate new function
	feSharedPtr<feMeshFunction> subf(new feSimple1dMeshFunction(f1->Mesh));
	subf->Clear();


	//fill new function
	feVector sublayer(f1->Mesh->NodesCount());
	for(int i=0; i<f1->LayersCount(); i++)
	{
		//calc du/dt
		for(int j=0; j<f1->Mesh->NodesCount(); j++)		
			sublayer[j] = f1->GetU((int)j, (int)i) - f2->GetU((int)j, (int)i);

		//add layer
		subf->AddLayer(f1->GetTime(i), &sublayer);
	}

	return subf;
}


//---------------------------------------------------------------------------------------------------------------------------------
feSharedPtr<feMeshFunction> feSimple1dMeshFunction::DiffDxOperation::operator() (feMeshFunction *iFunction)
{
	//cast
	feSimple1dMeshFunction *f = dynamic_cast<feSimple1dMeshFunction*>(iFunction);
	if(!f) return 0;


	//allocate new function
	feSharedPtr<feMeshFunction> df(new feSimple1dMeshFunction(f->Mesh));
	df->Clear();


	//fill new function
	feVector diflayer(f->Mesh->NodesCount());
	for(int i=0; i<f->LayersCount(); i++)
	{
		//calc du/dx
		for(int j=0; j<f->Mesh->NodesCount(); j++)		
			diflayer[j] = f->GetDuDx((int)j, (int)i);

		//add layer
		df->AddLayer(f->GetTime(i), &diflayer);
	}

	return df;
}


//---------------------------------------------------------------------------------------------------------------------------------
feSharedPtr<feMeshFunction> feSimple1dMeshFunction::IntDtOperation::operator() (feMeshFunction *iFunction)
{
	//cast
	feSimple1dMeshFunction *f = dynamic_cast<feSimple1dMeshFunction*>(iFunction);
	if(!f) return 0;


	//allocate new function
	feSharedPtr<feMeshFunction> intf(new feSimple1dMeshFunction(f->Mesh));
	intf->Clear();


	//fill new function
	feVector intlayer(f->Mesh->NodesCount());
	for(int i=0; i<f->LayersCount()-1; i++)
	{
		//calc du/dt
		for(int j=0; j<f->Mesh->NodesCount(); j++)		
			intlayer[j] += (f->GetU((int)j, (int)i) + f->GetU((int)j, (int)i+1)) * 0.5 * (f->GetTime(i+1) - f->GetTime(i));

		
	}

	//add layer
	intf->AddLayer(0, &intlayer);

	return intf;
}


//---------------------------------------------------------------------------------------------------------------------------------
void feSimple1dMeshFunction::AddLayer(double iT, feVector *iFunc)
{
	//check size
	if(iFunc->GetSize() != Mesh->NodesCount()) return;


/*
	//find place in vector and insert
	std::vector<Layer *>::iterator i;
	for(i = Layers.begin(); i != Layers.end(); ++i)
	{
		if((*i)->Time > iT)
		{
			--i;
			Layers.insert(i, new Layer(iT, iFunc));		
			return;
		}
	}	

*/

	//haven't found
	Layers[ActLayersCount] = new Layer(iT, iFunc);
	ActLayersCount++;
}


//---------------------------------------------------------------------------------------------------------------------------------
int feSimple1dMeshFunction::GetUindex(fePoint *iP, int iStartIndex)
{
	mfMesh *mesh = (mfMesh*)Mesh;

	//check bounds
	if(iP->X <= mesh->Points[0].X)
	{
		return 0;
	}
	if(iP->X > mesh->Points[mesh->ActPointsCount-1].X)
	{
		return -1;
	}


	//check bounds 
	if(iStartIndex >= mesh->ActPointsCount) iStartIndex = mesh->ActPointsCount/2;


	int nodes_count = mesh->ActPointsCount;

	//lookup local
	int ind1=iStartIndex, ind2=iStartIndex;
	int res = -1;

	for(; ind1>=iStartIndex-2; ind1--, ind2++)
	{
		if(ind1 >= 0 && mesh->Points[ind1].X <= iP->X && mesh->Points[ind1+1].X >= iP->X)
		{
			return ind1;
		}
		else if(ind2 < nodes_count-1 && mesh->Points[ind2].X <= iP->X && mesh->Points[ind2+1].X >= iP->X)
		{
			return ind2;
		}
	}


	//binary search
	ind1 = 0;
	ind2 = nodes_count-1;
	while(ind1 < ind2)
	{
		int i = (ind1+ind2)/2;
		if(mesh->Points[i].X < iP->X)
			ind1 = i+1;
		else
			ind2 = i;
	}

	if(mesh->Points[ind1].X <= iP->X && mesh->Points[ind1+1].X >= iP->X)
		return ind1;
	return ind2;

	return -1;
}


//---------------------------------------------------------------------------------------------------------------------------------
int	feSimple1dMeshFunction::GetTindex(double iTime, int iStartIndex)
{
	//check bounds
	if(iTime <= Layers[0]->Time)
	{
	//	feObserverManager::GetInstance().Notify("WARNING: out of range looking for t = %e\n", iTime);
		return 0;
	}
	if(iTime >= Layers[ActLayersCount-1]->Time)
	{
		//feObserverManager::GetInstance().Notify("WARNING: out of range looking for t = %e\n", iTime);
		return -1;
	}

	
	//check bounds 
	if(iStartIndex >= ActLayersCount) iStartIndex = ActLayersCount/2;


	//find
	int i = -1;
	int ind1=iStartIndex, ind2=iStartIndex;
	for(; ind1 >=iStartIndex-2; ind1--, ind2++)
	{
		if(ind1 >= 0 && Layers[ind1]->Time <= iTime && Layers[ind1+1]->Time >= iTime)
		{
			return ind1;
		}
		else if(ind2 < ActLayersCount-1 && Layers[ind2]->Time <= iTime && Layers[ind2+1]->Time >= iTime)
		{
			return ind2;
		}
	}


	//binary search
	ind1 = 0;
	ind2 = LayersCount()-1;
	while(ind1 < ind2)
	{
		int i = (ind1+ind2)/2;
		if( Layers[i]->Time < iTime)
			ind1 = i+1;
		else
			ind2 = i-1;
	}

	if(Layers[ind1]->Time <= iTime && Layers[ind1+1]->Time >= iTime)
		return ind1;
	return ind2;


	//done
	return -1;
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetU(fePoint *iP, double iTime, int *iLookUp)
{
	//id for reading cache
	static int cacheid = 0;

	//get index using cache
	int i;
	if(iLookUp != 0) 
	{
		if(*iLookUp > 0)
			i = *iLookUp;
		else 
			i = GetUindex(iP, *iLookUp);
	}
	else
		i = GetUindex(iP, CacheInUse->Get(cacheid));


	if(i<0)	return 0;
	

	//set new item
	CacheInUse->Set(cacheid, i);
	

	//find fe start index
	if(i == Mesh->NodesCount()-1)
		i --;

	if(iLookUp) *iLookUp = i;


	//get point over interpolation
	double x[2] = { Mesh->GetNode(i).X, Mesh->GetNode(i+1).X};
	double f[2] = { GetU(i, iTime), GetU(i+1, iTime)};
	return feInterpolation::Linear(x, f, iP->X);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetU(fePoint *iP, int iLayerNum, int *iLookUp)
{
	if(iLayerNum > ActLayersCount) return 0;

	//id for reading cache
	static int cacheid = 1;

	//get index using cache
	int i = iLookUp != 0 ? GetUindex(iP, *iLookUp) : GetUindex(iP, CacheInUse->Get(cacheid));
	if(i<0)	{ return 0; }	

	//set new item
	CacheInUse->Set(cacheid, i);


	//find fe start index
	if(i == Mesh->ActPointsCount-1)
		i--;
	if(i==0) i++;

	if(iLookUp) *iLookUp = i;


	//get point over interpolation
	double x[3] = { Mesh->Points[i-1].X, Mesh->Points[i].X, Mesh->Points[i+1].X};
	double f[3] = {Layers[iLayerNum]->Func.Data[i-1], Layers[iLayerNum]->Func.Data[i], Layers[iLayerNum]->Func.Data[i+1]};
	return feInterpolation::Quadratic(x, f, iP->X);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetU(int iNode, double iTime)
{
	if(ActLayersCount == 0) return 0;
	if(ActLayersCount < 3)
		return Layers[0]->Func.Data[iNode];


	//id for reading cache
	static int cacheid = 2;

	//get index using cache
	int i = GetTindex(iTime, CacheInUse->Get(cacheid));
	if(i<0)
		return 0;

	//set new cache item
	CacheInUse->Set(cacheid, i);


	//check bounds
	if(i == 0) i++;
	else if(i == ActLayersCount-1) i--;


	//get times
	double t1 = Layers[i-1]->Time;
	double t2 = Layers[i]->Time;
	double t3 = Layers[i+1]->Time;


	//get func over interpolation
	double x[3] = { t1, t2, t3 };
	double f[3] = {Layers[i-1]->Func.Data[iNode], Layers[i]->Func.Data[iNode], Layers[i+1]->Func.Data[iNode] };

	//done
	return feInterpolation::Quadratic(x, f, iTime);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetU(int iNode, int iLayerNum)
{
	if(iLayerNum >= ActLayersCount) iLayerNum = ActLayersCount - 1;
	if(iNode >= Layers[iLayerNum]->Func.GetSize()) return 0;
	return Layers[iLayerNum]->Func.Data[iNode];
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDt(fePoint *iP, double iTime)
{
	//id for reading cache
	static int cacheid = 3;

	//get index using cache
	int i = GetTindex(iTime, CacheInUse->Get(cacheid));
	if(i<0)
		return 0;

	//set new item
	CacheInUse->Set(cacheid, i);


	//check bounds
	if(i == 0) i++;
	else if(i == ActLayersCount-1) i--;


	//get ts, fs
	double ts[3] = { Layers[i-1]->Time, Layers[i]->Time, Layers[i+1]->Time };
	double fs[3] = {GetU(iP, i-1), GetU(iP, i), GetU(iP, i+1)};


	//done
	return feInterpolation::QuadDeriv(ts, fs, iTime);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDt(fePoint *iP, int iLayerNum)
{
	int tindex = 1;

	//check bounds
	if(iLayerNum == LayersCount()-1)
	{
		--iLayerNum;
		tindex = 2;
	}		
	else if(iLayerNum == 0)
	{
		++iLayerNum;
		tindex = 0;
	}


	//get ts, fs
	double ts[3] = { Layers[iLayerNum-1]->Time, Layers[iLayerNum]->Time, Layers[iLayerNum+1]->Time }; 
	double fs[3] = {GetU(iP, iLayerNum-1), GetU(iP, iLayerNum), GetU(iP, iLayerNum+1)};


	//done
	return feInterpolation::QuadDeriv(ts, fs, ts[tindex]);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDt(int iNode, double iTime)
{	
	//id for reading cache
	static int cacheid = 4;

	//get index using cache
	int i = GetTindex(iTime, CacheInUse->Get(cacheid));
	if(i<0)
		return 0;

	//set new item
	CacheInUse->Set(cacheid, i);


	//check bounds
	if(i == 0) i++;
	else if(i == ActLayersCount-1) i--;


	//get ts, fs
	double ts[3] = { Layers[i-1]->Time, Layers[i]->Time, Layers[i+1]->Time }; 
	double fs[3] = {GetU(iNode, i-1), GetU(iNode, i), GetU(iNode, i+1)};


	//done
	return feInterpolation::QuadDeriv(ts, fs, iTime);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDt(int iNode, int iLayerNum)
{
	//check bounds
	int tind = 1;
	if(iLayerNum == ActLayersCount-1)
	{
		--iLayerNum;
		tind = 2;
	}
	else if(iLayerNum == 0)
	{
		++iLayerNum;
		tind = 0;
	}


	//get ts, fs
	double ts[3] = { Layers[iLayerNum-1]->Time, Layers[iLayerNum]->Time, Layers[iLayerNum+1]->Time }; 
	double fs[3] = {GetU(iNode, iLayerNum-1), GetU(iNode, iLayerNum), GetU(iNode, iLayerNum+1)};


	//done
	return feInterpolation::QuadDeriv(ts, fs, ts[tind]);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDx(fePoint *iP, double iTime)
{
	//id for reading cache
	static int cacheid = 5;

	//get index using cache
	int i = GetUindex(iP, CacheInUse->Get(cacheid));
	if(i<0)
		return 0;

	//set new item
	CacheInUse->Set(cacheid, i);


	//check bounds
	if(i == 0) i++;
	else if(i == Mesh->NodesCount()-1) i--;


	//get ts, fs
	double xs[3] = {  Mesh->GetNode(i-1).X,  Mesh->GetNode(i).X,  Mesh->GetNode(i+1).X }; 
	double fs[3] = {  GetU(i-1, iTime), GetU(i, iTime), GetU(i+1, iTime) }; 


	//done
	return feInterpolation::QuadDeriv(xs, fs, iP->X);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDx(fePoint *iP, int iLayerNum)
{
	//id for reading cache
	static int cacheid = 6;

	//get index using cache
	int i = GetUindex(iP, CacheInUse->Get(cacheid));
	if(i<0)
		return 0;

	//set new item
	CacheInUse->Set(cacheid, i);


	//check bounds
	if(i == 0) i++;
	else if(i == Mesh->NodesCount()-1) i--;


	//get ts, fs
	double xs[3] = {  Mesh->GetNode(i-1).X,  Mesh->GetNode(i).X,  Mesh->GetNode(i+1).X }; 
	double fs[3] = {  GetU(i-1, iLayerNum), GetU(i, iLayerNum), GetU(i+1, iLayerNum) }; 


	//done
	return feInterpolation::QuadDeriv(xs, fs, iP->X);
	
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDx(int iNode, double iTime)
{	
	int xind = 1;

	//check bounds
	if(iNode == Mesh->NodesCount()-1)
	{
		--iNode;
		xind = 2;
	}
	else if(iNode == 0)
	{
		++iNode;
		xind = 0;
	}


	//get ts, fs
	double xs[3] = {  Mesh->GetNode(iNode-1).X,  Mesh->GetNode(iNode).X,  Mesh->GetNode(iNode+1).X }; 
	double fs[3] = {GetU(iNode-1, iTime), GetU(iNode, iTime), GetU(iNode+1, iTime)};


	//done
	return feInterpolation::QuadDeriv(xs, fs, xs[xind]);
}


//---------------------------------------------------------------------------------------------------------------------------------
double feSimple1dMeshFunction::GetDuDx(int iNode, int iLayerNum)
{
	//check bounds
	int xind = 1;
	if(iNode == Mesh->NodesCount()-1)
	{
		--iNode;
		xind = 2;
	}
	else if(iNode == 0)
	{
		++iNode;
		xind = 0;
	}


	//get xs, fs
	double xs[3] = {  Mesh->GetNode(iNode-1).X,  Mesh->GetNode(iNode).X,  Mesh->GetNode(iNode+1).X }; 
	double fs[3] = {  Layers[iLayerNum]->Func.Data[iNode-1], Layers[iLayerNum]->Func.Data[iNode], Layers[iLayerNum]->Func.Data[iNode+1] }; 


	//done
	return feInterpolation::QuadDeriv(xs, fs, xs[xind]);
}


//---------------------------------------------------------------------------------------------------------------------------------
bool feSimple1dMeshFunction::SaveToFile(char *iFileName, char *iMeshFileName)
{
	//check pointers
	if(!iFileName || !iMeshFileName) return false;


	//try to open

	FILE *fp = fopen(iFileName, "wb");
	if(!fp) return false;


	//save stride
	unsigned int size = sizeof(feReal);
	fwrite(&size, sizeof(size), 1, fp);


	Mesh->SaveToFile(iMeshFileName);

	//save some mesh info
	size = (unsigned int)Mesh->NodesCount();
	fwrite(&size, sizeof(size), 1, fp);
	size = (unsigned int)LayersCount();
	fwrite(&size, sizeof(size), 1, fp);


	//save function	
	for(unsigned int i=0; i<size; i++)
	{
		//time
		fwrite(&Layers[i]->Time, sizeof(Layers[i]->Time), 1, fp);

		//write vector
		fwrite(&Layers[i]->Func.Data[0], sizeof(feReal), GetLayer(i)->GetSize(), fp);
	}


	//done
	fclose(fp);
	return true;
}


//---------------------------------------------------------------------------------------------------------------------------------
bool feSimple1dMeshFunction::LoadFromFile(char *iFileName)
{
	//check pointers
	if(!iFileName) return false;


	//try to open
	std::ifstream fp(iFileName, std::ios::in | std::ios::binary);
	if(!fp.is_open()) return false;


	//read stride
	unsigned int stride = 0;
	fp.read((char*)&stride, sizeof(stride));
	
	//loading failed
	if(stride != sizeof(feReal)) return false;


	unsigned int ncount = 0;
	
	//load mesh info
	fp.read((char*)&ncount, sizeof(ncount));
	
	//loading failed
	if(ncount != Mesh->NodesCount()) return false;

	unsigned int lcount = 0;
	fp.read((char*)&lcount, sizeof(lcount));
	

	//resize vector
	ActLayersCount = lcount;


	//load function	
	for(unsigned int i=0; i<lcount; i++)
	{
		Layers[i] = new Layer();
		if(!Layers[i])
			return false;

		//time
		fp.read((char*)&Layers[i]->Time, sizeof(Layers[i]->Time));


		//read vector
		Layers[i]->Func.Allocate(ncount);
		fp.read((char*)&Layers[i]->Func.Data[0], sizeof(feReal)*ncount);
	}


	//done
	fp.close();
	return true;
}


//---------------------------------------------------------------------------------------------------------------------------------
void feSimple1dMeshFunction::Clear()
{
	for(unsigned int i=0; i<ActLayersCount; i++)
		delete Layers[i];
	ActLayersCount = 0;
}